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Abstract. In this paper, we present a new variational integrator for problems in Lagrangian mechanics. 
Using techniques from Galerkin variational integrators, we construct a scheme for numerical integration that 
converges geometrically, and is symplectic and momentum preserving. Furthermore, we prove that under 
appropriate assumptions, variational integrators constructed using Galerkin techniques will yield numerical 
methods that are in a certain sense optimal, converging at the same rate as the best possible approximation 
in a certain function space. We further prove that certain geometric invariants also converge at an optimal 
rate, and that the error associated with these geometric invariants is independent of the number of steps 
taken. We close with several numerical examples that demonstrate the predicted rates of convergence. 



There has been significant recent interest in the development of structure-preserving numerical methods 
for variational problems. One of the key points of interest is developing high-order symplectic integrators 
for Lagrangian systems. The generalized Galerkin framework has proven to be a powerful theoretical and 
practical tool for developing such methods. This paper presents a high-order Galerkin variational integrator 
that exhibits geometric convergence to the true flow of a Lagrangian system. In addition, this method is 
symplectic, momentum-preserving, and stable even for very large time steps. 

Galerkin variational integrators fall into the general framework of discrete mechanics. For a general 
and comprehensive introduction to the subject, the reader is referred to Marsden and West Q15]- Discrete 
mechanics develops mechanics from discrete variational principles, and, as Marsden and West demonstrated, 
gives rise to many discrete structures which are analogous to structures found in classical mechanics. By 
taking these structures into account, discrete mechanics suggests numerical methods which often exhibit 
excellent long term stability and qualitative behavior. Because of these qualities, much recent work has been 
done on developing numerical methods from the discrete mechanics viewpoint. See, for example, Hairer et al. 
[5] for a broad overview of the field of geometric numerical integration, and Marsden and West [TS]; Miiller 
and Ortiz |22| : Patrick and Cuell [23] discuss the error analysis of variational integrators. Various extensions 
have also been considered, including, Lall and West [TU]; Leok and Zhang [T7] for Hamiltonian systems; 
Fetecau et al. [8] for nonsmooth problems with collisions; Lew et al. [18]; Marsden et al. [20] for Lagrangian 
PDEs; Cortes and Martinez [5]; Fedorov and Zenkov [7]; McLachlan and Perlmutter [21 for nonholonomic 
systems; Bou-Rabee and Owhadi [5] [3J for stochastic Hamiltonian systems; Bou-Rabee and Marsden [T]; Lee 
et al. [TH [TJJ for problems on Lie groups and homogeneous spaces. 

The fundamental object in discrete mechanics is the discrete Lagrangian Ld : Q x Q x K — > K, where Q 
is a configuration manifold. The discrete Lagrangian is chosen to be an approximation to the action of a 
Lagrangian over the time step [0, h], 



or simply Ld (qo, q\) when h is assumed to be constant. Discrete mechanics is formulated by finding stationary 
points of a discrete action sum based on the sum of discrete Lagrangians, 



For Galerkin variational integrators specifically, the discrete Lagrangian is induced by constructing a discrete 
approximation of the action integral over the interval [0, h] based on a finite-dimensional function space and 
quadrature rule. Once this discrete action is constructed, the discrete Lagrangian can be recovered by solving 
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for stationary points of the discrete action subject to fixed endpoints, and then evaluating the discrete action 
at these stationary points, 

m 

(1) L d (q Q ,q 1 ,h) = ext h}bjL(q(cjh),q(cjh)). 

q„eU n {[0,h],Q) ~* 
g«(0)=?o,g„(/t)=«i J - 1 

Because the rate of convergence of the approximate flow to the true flow is related to how well the discrete 
Lagrangian approximates the true action, this type of construction gives a method for constructing and 
analyzing high-order methods. The hope is that the discrete Lagrangian inherits the accuracy of the function 
space used to construct it, much in the same way as standard finite-element methods. We will show that 
for certain Lagrangians, Galerkin constructions based on high-order approximation spaces do in fact result 
in correspondingly high order methods. 

Significant work has already been done constructing and analyzing these types of Galerkin variational 
integrators. In Leok [TJ], a number of different possible constructions based on the Galerkin framework are 
presented. In Leok and Shingel [15 , Hermite polynomials are used to construct globally smooth high-order 
methods. What separates this work from the work that precedes it is the use of a spectral approximation 
paradigm, which induces methods that exhibit geometric convergence. This type of convergence is established 
theoretically and demonstrated through numerical examples. 



1.1. Discrete Mechanics. Before discussing the construction and convergence of spectral variational inte- 
grators, it is useful to review some of the fundamental results from discrete mechanics that are used in our 
analysis. We have already introduced the discrete Lagrangian L d : Q x Q x R —> E, 

Ld (<?o,<?i, h) fa ext / L(q,q)dt. 

q ec 2 ([o,h],Q) J 
q(o)=qo,q(h)=qi 

and the discrete action sum, 

n-l „t 2 
§ ({%}Ll) = ^2 L d{qk,Qk+l) ~ / L(q,q)dt. 
k=l Jtl 

Taking variations of the discrete action sum and using discrete integration by parts leads to the discrete 
Euler-Lagrange equations, 

(2) D 2 L d (q k -i,qk) + DiL d {q k ,q k+ i) = 0, 

where D\ denotes differentiation with respect to the first argument and D 2 denotes differentiation with 
respect to the second argument. Given (q k -i,q k ), these equations implicitly define an update map, known 
as the discrete Lagrangian flow map, F^ d : Q x Q — >• Q x Q, given by Fi, d (qk-iiQk) = (<Zfc,?fc+i), where 
(qk-i,Qk) ! (qk,qk+i) satisfy Q. Furthermore, the discrete Lagrangian defines the discrete Legendre trans- 
forms, ¥ ± L d : Q x Q -> T*Q: 

¥ + L d : (qo,qi) -> (qi,Pi) = (qi,D 2 L d (q a ,qi)) , 
^~L d : (q ,qi) -> {qo,Po) = (qo, -DiLd (q ,qi)) ■ 

Using the discrete Legendre transforms, we define the discrete Hamiltonian flow map, F^ d : T*Q — > T*Q, 

F Ld ■ (qo,Po) -> Pi) = r+L d ((F-Ld) -1 (q , Po )) . 
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The following commutative diagram illustrates the relationship between the discrete Hamiltonian flow map, 
discrete Lagrangian flow map, and the discrete Legendre transforms, 




We now introduce the exact discrete Lagrangian , 

L d (q ,qi,h)= ext / L(q,q)dt. 

9(0)=9fc Jo 
l(h)=q k +i q qi 

q ec 2 ([a : h],Q) 

An important theoretical result for the error analysis of variational integrators is that the discrete Hamil- 
tonian and Lagrangian flow maps associated with the exact discrete Lagrangian produces an exact sampling 
of the true flow, as was shown in Marsden and West [19] . Using this result, Marsden and West [TO] shows 
that there is a fundamental relationship between how well a discrete Lagrangian Ld approximates the exact 
discrete Lagrangian and how well the corresponding discrete Hamiltonian flow maps, discrete Lagrangian 
flow maps and discrete Legendre transforms approximate each other. Since the exact discrete Lagrangian 
produces an exact sampling of the true flow, this in turn leads to the following theorem regarding the error 
analysis of variational integrators, also found in Marsden and West [T9"] : 

Theorem 1.1. (Variational Error Analysis) Given a regular Lagrangian L and corresponding Hamiltonian 
H, the following are equivalent for a discrete Lagrangian Ld-' 

(1) the discrete Hamiltonian flow map for Ld has error O (h p+1 j, 

(2) the discrete Legendre transforms of Ld have error O (h p+ J, 

(3) Ld approximates the exact discrete Lagrangian with error O (ft, p+1 ) . 

We will make extensive use of this theorem later when we analyze the convergence of spectral variational 
integrators. 

In addition, in Marsden and West |19j . it is shown that integrators constructed in this way, which are 
referred to as variational integrators, have significant geometric structure. Most importantly, variational 
integrators always conserve the canonical symplectic form, and a discrete Noether's Theorem guarantees 
that a discrete momentum map is conserved for any continuous symmetry of the discrete Lagrangian. The 
preservation of these discrete geometric structures underlie the excellent long term behavior of variational 
integrators. 

2. Construction 

2.1. Generalized Galerkin Variational Integrators. The construction of spectral variational integrators 
falls within the framework of generalized Galerkin variational integrators, discussed in Leok [T3] and Marsden 
and West [19]. The motivating idea is that we try to replace the generally non-computable exact discrete 
Lagrangian L^ (qk,qk+i) with a computable discrete analogue, L^f (qk,qk+i)- Galerkin variational integra- 
tors are constructed by using a finite-dimensional function space to discretize the action of a Lagrangian. 
Specifically, given a Lagrangian L : TQ —> R, to construct a Galerkin variational integrator: 

(1) choose an n-dimensional function space M n ([0, h] , Q) C C 2 ([0, h) , Q), with a finite set of basis 
functions (*)}™ =1 , 

(2) choose a quadrature rule Q (•) : F([0,/i],R) R, so that Q (/) = hJ2J=ibjf (cjh) J Q ' 1 / (t) dt, 
where F is some appropriate function space, 
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Figure 1. A visual schematic of the curve q n (t) € M" ([0, h],Q). The points marked with 
(x) represent the quadrature points, which may or may not be the same as interpolation 
points dih. In this figure we have chosen to depict a curve constructed from interpolating 
basis functions, but this is not necessary in general. 



and then construct the discrete action ({<$}•_■,) : Il™=i Qi ~ * ^ ( n °t to be confused with the discrete 
action sum S ({?fc}^U)), 

(/ n n \ \ m / n n \ 

l E (*) . E (*) = ft E ^ E fa* M) . E ' 
Vi=l i=l // j'=l \i=l f=l / 

where we use superscripts to index the weights associated with each basis function, as in Marsden and West 

ma 

Once the discrete action has been constructed, a discrete Lagrangian can be induced by finding stationary 
points q n (t) = Yli=i 1k4>i (*) °f the action under the conditions q n (0) = X)™=i (0) = Ik and q n (h) — 
Yh=i Qk&i ( h ) = for some given q k and q k+1 , 

m 

L d (q k ,q k+1 ,h) = ext § d ( {q\}™_, ) = cxt hS^bjL (q n (cjh) ,q n (cjh)) . 

9n(0)=9fc V / q n (0)=q k ~. 

q n (h)=q k + 1 qn(h)=q k+1 J 

A discrete Lagrangian flow map that result from this type of discrete Lagrangian is referred to as a Galerkin 
variational integrator. 



2.2. Spectral Variational Integrators. There are two defining features of spectral variational integrators. 
The first is the choice of function space M™ ([0, h] ,Q), and the second is that convergence is achieved not by 
shortening the time step h, but by increasing the dimension n of the function space. 



2.2.1. Choice of Function Space. Restricting our attention to the case where Q is a linear space, spectral 
variational integrators are constructed using the basis functions (pi (t) = Zj (t), where U (t) are Lagrange inter- 
polating polynomials based on the points hi — | cos (— ) + ^ which are the Chebyshev points ti — cos (^), 
rescaled and shifted from [—1, 1] to [0, h]. The resulting finite dimensional function space M" ([0, h] , Q) is 
simply the polynomials of degree at most n on Q. However, the choice of this particular set of basis functions 
offer several advantages over other possible bases for the polynomials: 

(1) the restriction on variations Yli=i (0) = Yl7=i ^ll^i W = reduces to 8q\ = 5q k = 0, 

(2) the condition q n (0) = q k reduces to q\ = q k , 

(3) the induced numerical methods have generally better stability properties because of the excellent 
approximation properties of the interpolation polynomials at the Chebyshev points. 
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Using this choice of basis functions, for any chosen quadrature rule, the discrete Lagrangian becomes, 

rn 

La{qk,qk+uh) = ext hS^bjL (q n (cjh) ,q n (cjh)) . 

Requiring the curve q n (t) to be a stationary point of the discretized action provides n — 2 internal stage 
conditions: 

m f dL ■ dL \ 

(3) h /] b i y-g^ (tin {cjh) , q n (cjh)) § v (cjh) + (fti (cjh) , q n (cjh)) <\> v {cjh)j = 0, p = 2, 

Combining these internal stage conditions with the discrete Euler-Lagrange equations, 

DiL d (qk-i,qk) + D 2 L d (q k ,q k+1 ) = 0, 
and the continuity condition q\. = q k yields the following set of n nonlinear equations, 

q\ = Ik, 

m f dL ■ dL \ 

foy^fy (q n (cjh),q n (cjh)) (j> p (cjh) + — (q n (cjh),q„ {c 3 h)) ij> p (cjh)j = 0, p = 2, ...,n- 1, 

m f dL - dL \ 

which must be solved at each time step k, and the momentum condition: 

m ( dL ■ dL \ 

hj^bj y— (q n (cjh) ,q n (cjh)) <p n {c h) + — (q n (cjh) ,q n {cjh)) <fin (cjh)J = p k , 

which defines Q for the next time step. Evaluating q k +i = q n (h) defines the next step for the discrete 
Lagrangian flow map: 

F Ld (ftb-i,fts) = (<7fc,<?fc+i) , 
and because of the choice of basis functions, this is simply q k +i = q k - 

2.2.2. n-Refinement. As is typical for spectral numerical methods (see, for example, Boyd [4]; Trefethen 
|25j). convergence for spectral variational integrators is achieved by increasing the dimension of the function 
space, M n ([0,h] ,Q). Furthermore, because the order of the discrete Lagrangian also depends on the order 
of the quadrature rule Q, we must also refine the quadrature rule as we refine n. Hence, for examining 
convergence, we must also consider the quadrature rule as a function of n, Q n . Because of the dependence 
on n instead of h, we will often examine the discrete Lagrangian L d as a function of Q x Q x N, 

TTln 

L d (q k , q k+u n) = Q) Sn (i (q n (*) , L (*))) = qn£U Sf^ h] Q) h J2 b ^ L (?" ( c ^ h ) ' ( c ^ h )) ' 

as opposed to the more conventional 

m 

L d (q k ,q k+ x,h) = q ^Jft MQ) Q (L (q n (t) , ft. (t))) = * £ ^ (ft. , ft. (c») . 



This type of refinement is the foundation for the exceptional convergence properties of spectral variational 
integrators. 
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3. Existence, Uniqueness and Convergence 

In this section, we will discuss the major important properties of Galerkin variational integrators and 
spectral variational integrators. The first will be the existence of unique solutions to the internal stage 
equations ([3| for certain types of Lagrangians. The second is the convergence of the one-step map that results 
from the Galerkin and spectral variational constructions, which will be shown to be optimal in a certain 
sense. The third and final is the convergence of continuous approximations to the Euler-Lagrange flow which 
can easily be constructed from Galerkin and spectral variational integrators, and the behavior of geometric 
invariants associated with the approximate continuous flow. We will show a number of different convergence 
results associated with these quantities, which demonstrate that Galerkin and spectral variational integrators 
can be used to compute continuous approximations to the exact solutions of the Euler-Lagrange equations 
which have excellent convergence and geometric behavior. 

3.1. Existence and Uniqueness. In general, demonstrating that there exists a unique solution to the 
internal stage equations for a spectral variational integrator is difficult, and depends on the properties of the 
Lagrangian. However, assuming a Lagrangian of the form 

L{q,q) = \q T Mq-V{q), 

it is possible to show the existence and uniqueness of the solutions to the implicit equations for the one-step 
method under appropriate assumptions. 

Theorem 3.1. (Existence and Uniqueness of Solutions to the Internal Stage Equations) Given a Lagrangian 
L:TQ^M. of the form 

L(q,q) = ^q T Mq-V(q), 

i/VV is Lipschitz continuous, bj > for every j and YHtLi^j = 1> an d M is symmetric positive-definite, 
then there exists an interval [0, h] where there exists a unique solution to the internal stage equations for a 
spectral variational integrator. 

Proof. We will consider only the case where q (t) £ K, but the argument generalizes easily to higher dimen- 
sions. To begin, we note that for a Lagrangian of the form, 

1 



L(q,q) = ^q T Mq-V(q) 



the equations 



Qk =Qk, 



m f dL ■ dL \ 

h)J)j y— (q n (cjh) ,q n (cjh)) (j) p {cjh) + — (q n (cjh) , <7„ (9/1)) P (cj7i)J =0, p = 2,...,n-l, 

m f dL ■ dL \ 

h Yl b 3 (*» ( c J h ) > <ln (Cjh)) 4>i (Cjh) + -gv (q~n (Cjh) , q n (cjh)) 4>i (cjh)J =Pk-l, 

take the form 

(5) Ax£ - f (q l ) = 0, 

where q l is the vector of internal weights, q % = {q\, q\, q r ^) T , A is a matrix with entries defined by 

(6) Ai,i =1, 

(7) A M =0, i = 2,...,n, 



(8) A p>i =h bjM<j>j (cjh) (j> p (cjh) , p = 2, ..., n;i = 1, n, 

3=1 
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and / is a vector valued function defined by 



/ (?) = 



Qk \ 

h T,7=l 6 3 W (Th=1 <fk<t>i ( c j h )) <f>n-l 

Pk-l I 



It is important to note that the entries of A depend on h. For now we will assume A is invertible, and that 
< ||j4^~ I] , for where Ai is the matrix A generated on the interval [0, 1]. Of course, the properties of 
A depend on the choice of basis functions {</>i}™ =1 , but we will establish these properties for the polynomial 
basis later. Defining the map: 

$ (V) = A~ l f (<f) ! 

it is easily seen that (JHj) is satisfied if and only if q l = $ [q % ), that is, q % is a fixed point of $(•)• If we 
establish that <& (•) is a contraction mapping, 

||$ (w l ) - $ (v l )\\ <k\\w l -v l \\ , 

II v / v/lloo — II lloo 

for some k < 1, we can establish the existence of a unique fixed point, and thus show that the steps of the 
one step method are well-defined. Here, and throughout this section, we use ||-|| to denote the vector or 
matrix p-norm, as appropriate. 

To show that <E> (•) is a contraction mapping, we consider arbitrary w l and v l : 

||$ («;*) - $ («*) L = II^T 1 / (w l ) - A- 1 / (v>) 1^ 

= ll^-M/M -/(«*)) L 
<II^1J|/K)-/K)L- 

Considering ||/ (u> 1 ) — / (w*) || , we see that 



(9) 11/ M "/(»*) L = 



3=1 



W £ (c,-fc) - VF ( £ «^ (c,-fc) 



for some appropriate index p* , Note that the first and last terms of ||/ (w l ) — f (v l ) || will vanish, so the 
maximum element must take the form of (Joj) . Let <\> % (t) — (<fri (t) ,<f>2(t),...,<fi n (t)). Let Cl be the Lipschitz 
constant for W(q). Now 



|/(O-/0/)| 



< 



3=1 

m 
3=1 



w ( £ ^ ^ c ^)J ^ w (£ ( c ^J 

/ n \ / n " 



|^p» (cj7i)| 



vi=l 



3=1 
m 

3=1 



|^p. (Cj/l)| 



£Ht-«£)& M) 



i=l 



|0p. (Cj/l)| 



< hJ2b,C L \\w* - ||^ (c^lL M)| 

3=1 
m 

max (||0 l (cj /i) 1 1 j \4> P ' (cjh)\) \\w l — 

3=1 J 
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= KCl max (||0* (cj/i)^ | V (cjh)\) \\w l - u 4 )^ . 
Hence, we derive the inequality 

11$ ( w *) - $ (« i )L < h ll^lloo ^ m f x (11^ M)lli M)D Ik* - vi \L > 

and since by assumption || 1 1| < ||.Aj~ || , 

||$ («,*) - $ (^)|| co < /i C L max (||^ (c^)^ |V (Cjh)\) ||«* " w*L • 

Thus if: 

^< (ll^rioo^maxdl^ (cj/i)^ |V M)|)) , 

then 

11$ (w l ) - $ (V)|| <fc||w l -w l || , 

II V / v/lloo — II lloo 

where k < 1, which establishes that $ (•) is a contraction mapping, and establishes the existence of a unique 
fixed point, and thus the existence of unique steps of the one step method. □ 

A critical assumption made during the proof of existence and uniqueness is that the matrix A is nonsin- 
gular. This property depends on the choice of basis functions 4>i. However, using a polynomial basis, like 
Lagrange interpolation polynomials, it can be shown that A is invertiblc. 

Lemma 3.1. (A is invertible) If {4>i}i = i is a polynomial basis of P n , the space of polynomials of degree at 
most n, M is symmetric positive-definite, and the quadrature rule is order at least 2n + 1, then A defined by 
^ - ([8| is invertible. 

Proof. We begin by considering the equation: 

Aq i = 0. 

Let q n (t) = Y^i=i Ik^i Considering the definition of A, Aq l = holds if and only if the following 
equations hold: 

9n (0) = 0, 

m 

(10) hY,bjMi l (c j h)^ p {c J h) = 0, p = 1, (n - 1). 

It can easily be seen that {<foj. is a basis of P n -i- Using the assumption that the quadrature rule is of 
order at least 2n — 1 and that M is symmetric positive-definite, we can see that (10) implies: 

/ Mi n (t) ^ (t) dt = 0, p = l,..,(n-l), 

Jo 

but, 

/ M?„ (t) fa (t) dt = 
Jo 

implies 

where (•, •) is the standard L 2 inner product on [0,h]. Since forms a basis for P n _i, q n G P n -ti 

and (•, •) is non-degenerate, this implies that q n (t) = 0. Thus, 

q n (0) = 

qn (t) = 



which implies that q n (t) = and hence q l = 0. Thus, Aq % = then q % = 0, from which it follows that A is 
non-singular. □ 



3.1 



bounded to prove Theorem 

We will do this by establishing 1 1 ^4_~ 1 1 1 oo — ll^i^lloo' wriere ^i i s ^ defined with h 



is 

oo 

is bounded. 



3.1 



Another subtle difficulty is that the matrix A is a function of h. Since we assumed that ||A _1 

we must show that for any choice of h, the quantity 1 1 ^4 1 

1. By Lemma 

we know that H^i^Hoo < 00 > which establishes the upper bound for ||j4 _1 || . This argument is easily 
generalized for a higher upper bound on h. 

Lemma 3.2. (II^ 1 )!^ < H^i^Hoo) For the matrix A defined by ^ - if h < 1, H^ -1 ^ < H^i^lL 
where A\ is A defined on the interval [0, 1] . 

Proof. We begin the proof by examining how A changes as a function of h. First, let {4>i}" = i be the basis 
for the interval [0, 1]. Then for the interval [0, h], the basis functions are 

t' 



and hence the derivatives are: 



(*) 



# (*) = -A 



Thus, if A\ is the matrix defined by Q — ffify on the interval [0, 1], then for the interval [0, h], 

1 




A = 



ft-f(n-l)x(n-l) 



where I nxn is the nxn identity matrix. This gives 



A~ 



a: 







(n-l)x(n-l) 



which gives 



-4; 



hi 





(n-l)x(n-l) 



< A: 



hh. 





-l)x(r. 



which proves the statement. 



□ 



3.2. Order Optimal and Geometric Convergence. To determine the rate of convergence for spectral 
variational integrators, we will utilize Theorem 1 1 . 1 1 and a simple extension of Theorem |1.1| 

Theorem 3.2. (Extension of Theorem |1.1| to Geometric Convergence) Given a regular Lagrangian L and 
corresponding Hamiltonian H, the following are equivalent for a discrete Lagrangian Ld (qo,qi,n): 

(1) there exists a positive constant K , where K < 1, such that the discrete Hamiltonian map for has 
error 0(K n ), 

(2) there exists a positive constant K, where K < 1, such that the discrete Legendre transforms of Ld 
have error O (K n ), 

(3) there exists a positive constant K , where K < 1, such that Ld is equivalent to a discrete Lagrangian 
with error O (K n ). 

This theorem provides a fundamental tool for the analysis of Galerkin variational methods. Its proof is 
almost identical to that of Theorem and can be found in the appendix. The critical result is that the 
order of the error of the discrete Hamiltonian flow map, from which we construct the discrete flow, has the 
same order as the discrete Lagrangian from which it is constructed. Thus, in order to determine the order 
of the error of the flow generated by spectral variational integrators, we need only determine how well the 
discrete Lagrangian approximates the exact discrete Lagrangian. This is a key result which greatly reduces 
the difficulty of the error analysis of Galerkin variational integrators. 

Naturally, the goal of constructing spectral variational integrators is constructing a variational method 
that has geometric convergence. To this end, it is essential to establish that Galerkin type integrators 
inherit the convergence properties of the spaces which are used to construct them. The order optimality 
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result is related to the problem of r-convergence (see, for example, Dal Maso [5]), as the Galerkin discrete 
Lagrangians are given by extremizers of an approximating sequence of variational problems, and the exact 
discrete Lagrangian is the cxtrcmizcr of the limiting variational problem. The r-convergence of variational 
integrators was studied in Miiller and Ortiz [52], and our approach involves a refinement of their analysis. We 
now state our results, which establish not only the geometric convergence of spectral variational integrators, 
but also order optimality of all Galerkin variational integrators under appropriate smoothness assumptions. 

Theorem 3.3. (Order Optimality of Galerkin Variational Integrators) Given an interval [0, h] and a La- 
grangian L : TQ — > R, let q be the exact solution to the Euler- Lagrange equations subject to the conditions 
q (0) = go ond q{h) — q^, and let q n be the stationary point of a Galerkin variational discrete action, i.e. if 
L% : Q x Q x R -> M, 

m 

Ld(qo,qh,h)= ext § rf ({qt}™-!) = ext h}bjL(q n (cjh),q n (cjh)), 



q n (0)=q ,q n {h)=q h 



q n (0)=q o ,q n (h)=q h 



3=1 



the 



If: 



q n = argmin h V] bjL (q n (cjh) , q n (cjh)) . 

q„m n ([0,h],Q) 3 = i 
q n (0)=q ,qn(h)=q h 



(1) there exists a constant Ca independent of h, such that, for each h, there exists a curve q n G 
M™ ([0,h],Q), such that, 



(q n (t),U*)) -(«(*) »?(<)) 



< C A h n , 



(2) there exists a closed and bounded neighborhood U C TQ, such that (q (t) , q (t)) € U , (q n (t) . q n (t) 
U for all t, and all partial derivatives of L are continuous on U , 

(3) for the quadrature rule Q (/) = h i bjf (cjh) w J Q f (t) dt, there exists a constant C g , such that, 



L {q n (t) , q n (t)) dt - h^bjL (q n (cjh) , q n (cjh)) 

3 = 1 

for anyq n eM n {[0,h},Q), 
(4) and the stationary points q, q n minimize their respective actions, 



< C g h 



n+l 



then 



I if (90, q h ,h)- L G d (g , q h ,h)\< C op h n +\ 



for some constant C op independent of h, i.e. discrete Lagrangian Ld has error O and hence the 

discrete Hamiltonian flow map has error O (/i n+1 ) . 

Proof. First, we rewrite both the exact discrete Lagrangian and the Galerkin discrete Lagrangian: 

r h 

L% (qo,qh,h) — L G (qo,qh,h)\ 



L(q(t),q(t))dt-g(L (g„ (*),«„ (t))) 
/ L (q (t ) , q (t ) ) dt - h V b 3 L (q n (c 3 h) , q n (c 3 h) ) 

/ L (q, q) dt - h V] bjL (q n , q r 
Jo „._i 



3 = 1 



10 



where in the last line, we have suppressed the t argument, a convention we will continue throughout the 
proof. Now we introduce the action evaluated on the q n curve, which is an approximation with error O (h n ) 
to the exact solution q: 



rh m ph rt, 

\ L (q, q) dt - h V] b 3 L (q n , q n ) = / L{q,q)dt~ / 
Jo - =1 Jo Jo 



(11a) 



(lib) 



< 



(qnA-n^j dt 

/ L f q n , q n J dt-h^2 b 3 L (<ln, <?« 

L q)dl - j L (<7„,<?„) ( 

/ L (q n ,q n ) dt - hS^bjL (q n ,q n ) 
Jo v ' j=1 



Considering the first term (11a): 

rh 



o 



L(q,q)dt L(q n ,q n )dt 



< 



(q,q) - L (q n ,qn) 



dt 



dt. 



By assumption, all partials of L are continuous on U, and since U is closed and bounded, this implies L 
is Lipschitz on U. Let L a denote that Lipschitz constant. Since, again by assumption, (q, q) G U and 

<ln, q-n ) € U, we can rewrite: 



L(q,q) ~ L (q n ,q, 



dt< L a 



(q,q) - (?n,<?™) 



dt 



< / L a C A h n dt 
Jo 

= L a C Ah n+1 , 

where we have made use of the best approximation estimate. Hence, 

< L a dh n+1 . 



(12) 



L(q,q)dt- / L(q n ,q n ) dt 



Next, considering the second term (lib), 



/ L[q n ,q n \ dt-h22 b 3 L (^n . £ 

J ■ i 



since g„, the stationary point of the discrete action, minimizes its action and q n £ M" ([0, h] , Q), 

mm./, 
(13) hJ2 b 3 L {qn, L) < h^bjL < / L [q n , 4n) dt + C g h n+1 

3 = 1 j=l J ° 

where the inequalities follow from the assumptions on the order of the quadrature rule. Furthermore, 

l-h 



h V bjL (q n ,hn)> L (q n , £„) dt - C g h n+1 

3 = 1 J ° 

> / L(q,$)dt-C g h n+1 
Jo 
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(14) 



> / L(q n ,q- n )dt-L a C A h n+1 -C g h 



n+l 



where the inequalities follow from ( 12 1, the order of the quadrature rule, and the assumption that q minimizes 
its action. Putting ( 13 ) and ( 14 1 together, we can conclude: 

(15) 



if \ 

/ L (q n , q n ) dt - h}^ bjL (q 

ni qn 



< {L a C A + Cg)h 



n+l 



Now, combining the bounds ( |12[ ) and (151 in (11a) and (lib), we can conclude 

\Lf (q ,q h ,h)-L% (q ,q h ,h)\ < (2L a C A + C g ) h n+1 
which, combined with Theorem |1.1[ establishes the order of the error of the integrator. 



□ 



The above proof establishes a significant convergence result for Galerkin variational integrators, namely 
that for sufficiently well behaved Lagrangians, Galerkin variational integrators will produce discrete ap- 
proximate flows that converge to the exact flow as h — > with the highest possible order allowed by the 
approximation space, provided the quadrature rule is of sufficiently high order. 

We will discuss assumption 4 in j ]3.3| While in general we cannot assume that stationary points of the 
action are minimizers, it can be shown that for Lagrangians of the canonical form 

L(q,q) = q T Mq-V(q), 

under some mild assumptions on the derivatives of V and the accuracy of the quadrature rule, there always 
exists an interval [0, h] over which stationary points are minimizers. In [3.3 we will show the result extends 
to the discretized action of Galerkin variational integrators. A similar result was established in Miiller and 
Ortiz [22] . 

Geometric convergence of spectral variational integrators is not strictly covered under the proof of order 
optimality. While the above theorem establishes convergence of Galerkin variational integrators by shrinking 
h, the interval length of each discrete Lagrangian, spectral variational integrators achieve convergence by 
holding the interval length of each discrete Lagrangian constant and increasing the dimension of the ap- 
proximation space M™ ([0, h] ,Q). Thus, for spectral variational integrators, we have the following analogous 
convergence theorem: 

Theorem 3.4. (Geometric Convergence of Spectral Variational Integrators) Given an interval [0, h] and a 
Lagrangian L : TQ — > R, let q be the exact solution to the Euler- Lagrange equations, and q n be the stationary 
point of the spectral variational discrete action: 



Lf(q ,q h ,n)= ^jf^ MM?=i) 

q n (0)=q ,q n (h)=q h 



ext h'S^b n L(q n {c n h),q n (cn 1 h)) 

9„£M"([0,/i].Q) J v 



If: 



(1) there exists constants C a ,Ka, Ka < 1 ; independent of n such that, for each n, there exists a curve 
q n € M" ([0,h],Q), such that, 

(q,q) - (q n ,L) < C A K% 

(2) there exists a closed and bounded neighborhood U C TQ, such that, (q (t) , q (t)) € U and (jj n (t) , q n (t) 
U for all t and n, and all partial derivatives of L are continuous on U , 

(3) for the sequence of quadrature rules Q n (f) — ^Zfci bnjf {c nj h) ~ J Q f (t) dt, there exists constants 
C g , Kg, K g < 1, independent of n such that 



L (q n (t) ,q n (t)) dt - h^b nj L (q n (c n] h) ,q n (c nj h)) 

3 = 1 

for any q n G M n ([0, h] , Q), 
(4) and the stationary points q, q n minimize their respective actions, 



< CgK-. 



then 

(16) \Lf (q , qi )-L s d (q ,qi,n)\ <C S K^ 

for some constants C S ,K S , K s < 1, independent of n, and hence the discrete Hamiltonian flow map has 
error O (K£). 

The proof of the above theorem is very similar to that of order optimality, and would be tedious to repeat 
here. It can be found in the appendix. The main differences between the proofs are the assumption of 
the sequence of converging functions in the increasingly high-dimensional approximation spaces, and the 
assumption of a sequence of increasingly high-order quadrature rules. These assumptions are used in the 
obvious way in the modified proof. 



3.3. Minimization of the Action. One of the major assumptions made in the convergence theorems (3.3) 



and (3.4) is that the the stationary points of both the continuous and discrete actions are minimizers over 
the interval [0, ft]. This type of minimization requirement is similar to the one made in the paper on T- 
convergence of variational integrators by Miiller and Ortiz 22 . In fact, the results in Muller and Ortiz [52] 
can easily be extended to demonstrate that for sufficiently well-behaved Lagrangians of the form 

L{q,q) = \q T Mq~V{q), 

where q 6 C 2 ([0, h] ,Q), there exists an interval [0, h], such that stationary points of the Galerkin action are 
minimizers. 



Theorem 3.5. Consider a Lagrangian of the form 

L{q,q)= l -q T Mq-V(q) 

where q £ C 2 ([0, h] , Q) and each component q d of q, q d G C 2 ([0, h] , Q), is a polynomial of degree at most 
s. Assume M is symmetric positive-definite and all second-order partial derivatives of V exist, and are 
continuous and bounded. Then, there exists a time interval [0, h] such that stationary points of the discrete 
action, 

m , 

on this time interval are minimizers if the quadrature rule used to construct the discrete action is of order 
at least 2s + 1. 

We quickly note that the assumption that each component of q, q d , is a polynomial of degree at most s 
allows for discretizations where different components of the configuration space are discretized with polyno- 
mials of different degrees. This allows for more efficient discretizations where slower evolving components 
are discretized with lower-degree polynomials than faster evolving ones. 

Proof. Let q n be a stationary point of the discrete action (•), and let Sq be an arbitrary perturbation of 
the stationary point q n , under the conditions Sq d G Ps d i <5q (0) = Sq(h) — 0, which is uniquely defined by 
C Q. Then, 



^{{ql + 5ql) n i=1 )-^ d [{ql} n i=1 

m ( 1 \ ™ / 1 

m / 1 i 

= h J2 b A o + 6< i) T M (?» + 5( *) - V (?» + W - oQnMq~ n + V (q n ) 



J 



Making use of Taylor's remainder theorem, we expand: 



V (q n + Sq) = V {q n ) + W (q n ) ■ Sq + ^SqlRSq n , 
13 



where \R im \ < supj . 



d 2 v 



. Using this expansion, we rewrite 



®d (U + hlY %=1 ) - §d (Ul}tl) = h J2 b > \ 2 & + 5 ^f M & + 5q)-V (q n ) - W (q n ) ■ Sq 



1 • 



- -Sq 1 RSq - i^-qi Mq n + V (q n ) 
which, given the symmetry in M, rearranges to: 

m / 1 1 \ 

3 ' 

Now, it should be noted that the stationarity condition for the discrete Euler-Lagrange equations is 

m 

h ]T b 3 (%M6q - W (q n ) -5q)=0 
for arbitrary Sq, which allows us to simplify the expression to 

^ (H + KYi,) ^ ({«Ur=i) = 6 i (\&fM8q \8q T R8q 



Now, using the assumption that the partial derivatives of V are bounded, \Ri m \ < < C R , an d 

standard matrix inequalities, we get the inequality: 

(17) 5g T iMg < ||itfg|| 2 \\6q\\ 2 < ||i?|| 2 INI* < \\R\\ F \\6qf 2 < DC R \\5qf 2 = DC R Sq T Sq, 

where D is the number of spatial dimensions of Q. Thus 

m , \ m / 1 1 \ 

ft I> ( 2*9 TM<y 9- 2*« TiW «) ^ h Y< hj U 5qTM5q ~2 DCR5qT5q ) ■ 

3 3 

Because M is symmetric positive-definite, there exists m > such that x T Mx > mx T x for any x. Hence, 
( 2^ T M<S9- -DCrfM > /»X)6j f -mSq T Sq- -DC R 6q T Sq\ . 

3 3 

Now, we note that since each component of 5q is a polynomial of degree at most s, 6q T Sq and Sq T Sq are both 
polynomials of degree less than or equal to 2s. Since our quadrature rule is of order 2s + 1, the quadrature 
rule is exact, and we can rewrite 

m , \ \ ( h 

h^bj \^-m5q T 5q - -DC R 5q T Sqj = -J mSq T 5q - DC R Sq T Sqdt 

= \ m5 f 5 i At - J q DC R Sq T Sqdtj . 

From here, we note that Sq G i?d ([0, h] ,Q), and make use of the Poincare inequality to conclude 
\{f mSq T 5qdt-J^ nC R Sq T Sqdtj > ^ (m^ J Sq T Sqdt — DC R J Sq T 5qdt\ 



U ^ -DC R )j\q-8qdt. 



2 V h 2 



Since J Q 5q T Sqdt > 0, 



s d ({«i+^}r=i)-s d ({«i}r=i 



> 



1 / mir 2 



2 
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DC R / Sq T Sqdt > 



so long as h< J . □ 
3.4. Convergence of Galerkin Curves and Noether Quantities. 

3.4.1. Galerkin Curves. In order to construct the one-step method, spectral variational integrators determine 
a curve, 

n 



»=i 



which satisfies 



q n (t) = argmin h^b 3 L (q n (cjh) , q n (cjh)) 



,6M"([0,k],Q) 

»(0)=gs,,gn(/i)=9fe+i 



3 = 1 



Evaluating this curve at h defines the next step of the one-step method, qk+i — q n (h), but the curve itself 
has many desirable properties which makes it a good continuous approximation to the true solution of the 
Euler Lagrange equations q(t). In this section, we will examine some of the favorable properties of q n (t), 
hereafter referred to as the Galerkin curve. 

However, before discussing the properties of the Galerkin curve, it is useful review the different curves 
with which we are working. We have already defined the Galerkin curve, q n (t), and we will also be making 
use of the local solution to the Euler-Lagrange equations q (t), where 

rh 



q(t)= argmin / L (q (t) , q (t)) dt. 
n pc 2 ao.h].o\ Jo 



q ec 2 ([o,h],Q) 

q(0)=q k ,q(h)=q k + 1 

However, while for each interval q satisfies the Euler-Lagrange equations exactly, it is not the exact solution 
of the Euler-Lagrange equations globally, as q^ ^ §kh (<7cb Qo)i where 3>j ((fa, 3o) is the flow of the Euler- 
Lagrange vector field. This is particularly important when discussing invariants, where the invariants of q 
remain constant within a time-step, but not from time-step to time-step. 

The first property of the Galerkin curve that we will examine is its rate of convergence to the true flow 
of the Euler-Lagrange vector field. There are two general sources of error that affect the convergence of 
these curves, the first being the accuracy to which the curves approximate the local solution to the Euler- 
Lagrange equations over the interval [0, h] with the boundary (q^, qk+i), and the second being the accuracy 
of the boundary conditions (ft+i) as approximations to a true sampling of the exact flow. Numerical 
experiments will show that often the second source of error dominates the first, causing the Galerkin curves 
to converge at the same rate as the one-step map. However, the accuracy to which the Galerkin curves 
approximate the true minimizers independent of the error of the boundary can also be established under 
appropriate assumptions about the action. Two theorems which establish this convergence are presented 
below. 

Before we state the theorems, we quickly recall the definitions of the Sobolev Norm 



iw^ffltyi]) _ I iI/Hlp([o,/i]) + 



LP([0,h]) 




We will make extensive use of this norm when examining convergence of Galerkin curves. 

Theorem 3.6. (Geometric Convergence of Galerkin Curves with n-Refinement) Under the same assumptions 
as Theorem \3.4\ if at q, the action is twice Frechet differentiable, and if the second Frechet derivative of the 
action D 2 & (•) [•, •] is coercive in a neighborhood U of q, that is, 

D 2 &{v) [6q,Sq]>C f \\dq\\ 

for all curves Sq £ Hq ([0, h] , Q) and all v G U , then the curves which minimize the discrete action converge 
to the true solution geometrically with n-refinement with respect to IHIwMffo />])• Specifically, if the discrete 
Hamiltonian flow map has error O (K™), K s < I, then the Galerkin curves have error O t\J K s ) . 
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Proof. We start with the bound (16) given at the end of Theorem 3.4 

\L% {dh,Qk+i)-L s d (gfc,gfc+i,n)| < C S K^ 

and expand using the definitions of (qk,Qk+i) and (q k , qk+i, n), as well as the assumed accuracy of 
the quadrature rule Q n to derive 

C a K n s > |Lf (q k ,q k+1 )-L% {q k ,q k+ i,n)\ 



(18) 



L(q,q)dt- h^b nj L (q n (c nj h) , q n (c nj h) ) 

3=1 



(19) 



which implies: 



> 



L (q, q) dt 



L(q n ,q n )<it 



C g K n g 



\&{q n )-&{q)\-C g K™ 

(C s + Cg) K n s > I© (q n ) - 6 (g)| 



because K s > (see the proof of Theorem 3.4 in the appendix). Using this inequality, we make use of a 
Taylor expansion of © (qn), 



&(q n )=&(q)+T)&(q) [q n ~q] 

for some v £ U, to see that 

(C s + C g )K2>\e(q n )-e(q)\ 

= © (q) + D© (q) [<7„ — q] - 



1 



D 2 6 (v) [q n -q,q n - q] 



1 



But 



DG (q) [q n - q] 



h dL 
dq 
h 'dL 



D 2 © (q) % 



dL 



q,q n 



6(q) 



{q, i) {q n - q) + -^r (q, q) [q n - i) dt 

q)dt 



M Voq { "-' l)] 



0, 



because q n (0) = q (0) and q n (h) — q (h) by definition (note that this implies (q n — q) £ Hq ([0, h] , Q)). Then 



(C s + C g )Kl l > |D 2 ©H % 
> C 



9l I 



/ n 2 

CyK s >\\qn — q\\ w i,i([Q.h}) 



where C 



C S +Cg 



□ 



This result shows that Galerkin curves converge to the true solution geometrically with n-refinement, 
albeit with a larger geometric constant, and hence a slower rate. By simply replacing the bounds (18 1 and 
(19) from Theorem 3.4 with those from Theorem 3.3 and the term C S K™ with C op h p , an identical argument 
shows that Galerkin curves converge at half the optimal rate with /i-refinement. 

Theorem 3.7. (Convergence of Galerkin Curves with h- Refinement) Under the same assumptions as The- 



orem 



3.3 if at q, the action is twice Frechet differentiable, and if the second Frechet derivative of the 



action D S (•) [•, •] is coercive with a constant Cf independent of h in a neighborhood U of q, for all curves 
Sq £ Hq ([0, h] , Q), then if the discrete Lagrange map has error O the Galerkin curves have error at 

* n ll"lliyi.i([o h])- IfCf is a function of h, this bound becomes O i^C f (h) '/stJ. 



most O ( ft. P 2 
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Like the requirement that the stationary points of the actions are minimizers, the requirement that the 
second Frechet derivative of the action is coercive may appear quite strong at first. Again, the coercivity will 
depend on the properties of the Lagrangian L, but we can establish that for Lagrangians of the canonical 
form, 

L{q,q) = \q T Mq-V{q), 
there exists a time step [0, h] over which the action is coercive on Hq ([0, h] ,Q). 
Theorem 3.8. (Coercivity of the Action) For Lagrangian of the form 

L(q,q) = ^q T Mq-V(q), 

where M is symmetric positive-definite, and the second derivatives of V (q) are bounded, there exists an 
interval [0, h] over which the action is coercive over Hq ([0, h] ,Q), that is, 

D 2 <S{v) [6q,6q]>C f \\Sq\\ 

for any Sq e Hq ([0, h] , Q) and any v £ C 2 ([0, h],Q). 

Proof. First, we note that if 

&iy) = J \v T Mv-V{v), 

then 

r h 

D 2 6 (v) [Sq, Sq] = / Sq T M8q - 5q T H (u) 6qdt 
Jo 

nh rh 

= / 5q T M5qdt- / 5q T H{v)5q&t 
Jo Jo 

where H (v) is the Hessian of V (is) at the point v. Since M is symmetric positive-definite, and the second 
derivatives of V (•) are bounded, then there exists C r and m such that: 



(20) 



Sq T MSqdt> / mSq T Sqdt 
Jo 

/ 5q T H{v)8qdt< I DC r Sq T Sqdt, 
Jo Jo 



(see (17) for a derivation of (20)). Hence, 

D 2 6 (y) [Sq, Sq] > / mSq T 8qdt - / DC r f T fdt 
Jo Jo 

1 fh i rh rh 

(21) = -m / Sq T Sqdt + ~m 5q T Sqdt - DC r / Sq T 6qdt. 

2 Jo 2 j Jo 

Considering the last two terms in (21), and noting that 8q € Hq ([0,h],Q), we make use of the Poincare 
inequality to derive: 

1 C h [ h m-K 2 f h r h 

-ml Sq T Sqdt-DC r 5q T 5qdt>—— Sq T Sqdt-nC r 6q T 5qdt 

2 Jo Jo 2ft J Jq 

(22) > f^-DCr) [ h Sq T Sqdt. 



Thus, substituting (p2| in for the last two terms of (21 ), we conclude: 



D 2 6 (,, 4) [Jj, 5q] > (j£ - DC T \ J Sg T Sgdt + | / 5q T 5q4t 
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Figure 2. Conserved and approximately conserved Noether quantities and the resulting 
constrained solution space. Suppose that both p T q = 1 and p 2 + q 2 = 5 were conserved 
quantities for a certain Lagrangian. Then the solutions of the Euler-Lagrange equations 
would be constrained to the intersections of these two constant surfaces in phase space; in 
the above diagram, this is the intersection of the dashed and solid lines. If these quantities 
were conserved up to a fixed error along a numerical solution, then the numerical solution 
would be constrained to the intersection of the shaded regions in the above figure. The 
constraint of the numerical solution to these regions is what leads to the many excellent 
qualities of variational integrators. 



DC, 



vj J (ini^qo^) + iNii2 ([(Vl]) ) , 



2 ' V 2h 2 

and making use of Holder's inequality, we see that ||<5g|| L 2([ h ^ > hi \\Sq\\ L i^ thus 



m j 


/ U17T 2 


T'l 


K 2h 2 


mh 


( m-K 2 


~2~' 


\ 2h 


mh 


f mir 2 


T' 


( 4h 



DC, 



> min — , — hDC, 



h\\Sq\\l W]) ) 

If y 
2 (H*«Hii([o,fc]) + H^HiVIo.fcljJ 



- hDC r j j ||<yg||^i.i([ ,fc]) 

which establishes the coercivity result. □ 



3.4.2. Noether Quantities. One of the great advantages of using variational integrators for problems in geo- 
metric mechanics is that by construction they have a rich geometric structure which helps lead to excellent 
long term and qualitative behavior. An important geometric feature of variational integrators is the preser- 
vation of discrete Noether quantities, which are invariants that are derived from symmetries of the action. 
These are analogous to the more familiar Noether quantities of geometric mechanics in the continuous case. 
We quickly recall Noether's theorem in both the discrete and continuous case, which will also help define the 
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notation used throughout the proofs that follow. The proofs of both these theorems can be found in Hairer 
et al. 0- 



Theorem 3.9. (Noether's Theorem) Consider a system with Hamiltonian H (p, q) and Lagrangian L (q, q). 
Suppose {g s : s £ R} is a one-parameter group of transformations which leaves the Lagrangian invariant. Let 

d 



as 



9s (?) 



be defined as the vector field with flow g s {q), referred to as the infinitesimal generator, and define the 
canonical momentum 

dL 

P= (ff,ff). 

Then 

i (p, q) = P T a (?) 

is a first integral of the Hamiltonian system. 

Theorem 3.10. (Discrete Noether's Theorem) Suppose the one-parameter group of transformations leaves 
the discrete Lagrangian Ld (qk,qk+i) invariant for all (qk,qk+i)- Then: 

pl +1 a(q k+1 ) =pla{q k ) 

where 

Pk = -L>iL d (q k ,q k+ i) , 
Pk+i = D 2 L d (q k ,q k+ i) ■ 

For the remainder of this section, we will refer to I(q,p) as the Noether quantity and p^a(q n ) = 
Pn+i a (qn+i) as the discrete Noether quantity. 

For Galerkin variational integrators, it is possible to bound the error of the Noether quantities along the 
Galerkin curve from the behavior of the analogous discrete Noether quantities of the discrete problem and, 
more importantly, this bound is independent of the number of time steps that are taken in the numerical 
integration. This is significant because it offers insight into the excellent behavior of spectral variational 
integrators even over long periods of integration. 

The proof of convergence and near preservation of Noether quantities is broken into three major parts. 
First, we note that on step k of a numerical integration the discrete Noether quantity arises from a function 
of the Galerkin curve and the initial point of the one-step map (qk-i,qk), and that a bound exists for the 
difference of this discrete Noether quantity evaluated on the Galerkin curve and evaluated on the local exact 
solution to the Euler-Lagrange equations q. Second, we show that a bound exists for the difference of the 
discrete Noether quantity on the local exact solution of the Euler-Lagrange equations and the value of the 
Noether quantity of the local exact solution, which is conserved along the flow of the Euler-Lagrange vector 
field. Finally, we show that under certain smoothness conditions, there exists a point-wise bound between 
the Noether quantity evaluated on the Galerkin curve and the Noether quantity evaluated on the local exact 
solution. Thus, we establish a point-wise bound between the Noether quantity evaluated on the Galerkin 
curve and the discrete Noether quantity, and a bound between the discrete Noether quantity and the Noether 
quantity, which leads to a point-wise bound between the Noether quantity evaluated on the Galerkin curve, 
and the Noether quantity which is conserved along the global flow of the Euler-Lagrange vector field. 

Throughout this section we will make the simplifying assumptions that 



q n = 2_^q l k$ 



where q\ ~ qk, and thus 

dq n _ 
dqk 

This assumption significantly simplifies the analysis. 
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We begin by bounding the discrete Noether quantity by a function of the local exact solution of the 
Euler-Lagrange equations. 

Lemma 3.3. (Bound on Discrete Noether Quantity) Define the Galerkin Noether map as: 

T 

'dL . i \ i dL 



— (q,q)4>i + (q,q)4>i 



Id (q(t) ,qk) = - ( h^bj 

3 = 1 

and note that the discrete Noether quantity is given by 

Id (tin, 4k) =Pn a (<lk 



a(q k 



Assuming the quadrature accuracy of Theorem (3.4) with n-refinement and Theorem (3.3) with h-refinement, 
if^(q,q), §^ (q, q) and o,re Lipschitz continuous, \\4>i\\ l°° ([o h]) * s bounded with n refinement, and 

llSn — 3llvfi.i([o h]) * s bounded below by the quadrature error, then 

\Id{q n ,qk) - Id(q,qk)\ < C\a{q k )\ (\\q n - q\\ w i,i([o, h ]) + Mn - s|lz,°°([o,fc]) + \% ~ $\\ L °°( [0ih]) 
for some C independent of n and h. 

Proof. We begin by expanding the definitions of the discrete Noether quantity: 

T 



\ld(q n ,qk) - id(q,qk)\ 



M5> 

3=1 



dL ._ ■_ , dL , i . • ' 

-Qq (q n ,q n ) 91 + -gv [q n ,gn) fa 



a{qk) 



3=1 



91 ,_ , dL , j.. • 

-q^ \q,q) fa + -qt [q,q) fa 



a(q k ) 



9g 



<9L 



dL 



dL 



dq 



dq 



< 



h lL, b i (?n,?n) - ¥ (?,?)J ^i - (q n ,L) ~ S (q,q)j fa 

4=1 



a (<j fc ) 
a{qk)\ ■ 



Now we introduce the function e q (•, •) which gives the error of the quadrature rule, and thus 



\Id(q n ,qk) ~ Id(q,qk)\ < 



h (dL, i 



dL 



(9, g) 



<9g 

e g (q n - q, q n - q 



dL - i , 91 , . 1 ■ 
gvC*.i9n)-gv(g,?)J0ld* 



k(gfe)l 



Integrating by parts, we get: 

- ^ (g,gfe)l < 



+ ^(gn,gn)-^(g,g)J 



>i&t 

a{qk)\ ■ 



Introducing the Lipschitz constants L\ for L2 for and L3 for 



l^(<7n,9fe) -id(g,gfe)| < ( / (Li +L3) \(q n ,q n ) - (g,g)| \fa\dt + 2L 2 [\\fa\\ L oo^ 0xh] ^ 

(\\q~n - g|lz,=o([0,fe]) + Mn ~ gllz^cfc])) + e ? (g*n - 9,tin ~ $) J |« (<?fc)l 
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< (Li + L 3 ) ||^i||ioo ([0)A]) |o(g*)l yj o \(qn,q n ) - (q,q)\dtj 

+ 2^2 (ll&llz»([o,fc])) l a («*)l (\Wn ~ q\\ L °°([a,h]) + \\9n ~ ?|| L oo ([0>h]) 

+ e q (q n ~q,q n - q) \a(q k )\ . 

We now make the simplification that the quadrature error \e q (•, -)| serves as a lower bound for \\q n — q\\w^A([o,h])- 
While this may not strictly hold, all of our estimates on the convergence for q n imply this bound, and hence 
it is a reasonable simplification for establishing convergence in this case. Now, note that ||<^i||j,ooqo fel) ^ s 
invariant under h rcscaling, and let 



C = max(L 1 + L 3) 2L 2 )\\cf> 1 \\ Lo 



i[o,h]) 



1 



to get 



\Id(Qn,qk)-Id(q,qk)\ < C\a(q k )\ [\\q„ - q\\ w ^([o,h]) + ^ n ~ 5lli«([o,h]) + \\$n ~ q\\ L ~ i[0ih]) 
which establishes the result. 



□ 



Lemma 3.3 establishes a bound between the discrete Noether quantity and Id (q,q k )- The next step is to 
establish a bound between 1^ (g, q k ) and the Noether quantity. 

Lemma 3.4. (Error Between Discrete Noether Quantity and True Noether Quantity) Assume that tfii (0) = 1 
and 4>i (h) = 0, and that the sequence {\a (qk)\}^—i * s bounded by a constant C a which is independent of N. 
Let 

P(t) = -gv(q(t),q(t))- 
Once again, let the error of the quadrature rule be given by e q (•,•). Then 



\l d (q,q k )-I(p(t),q(t))\<C a \e q (q,q)\ 



for any t £ [0, h]. 



Proof. First, we note that since q solves the Euler-Lagrange equations exactly, I (p(t) ,q{t)) is a conserved 
quantity along the flow, so it suffices to show the inequality holds for t = 0. We begin by expanding: 

h ( £ b ^ & & + ^ & 5) k ) « (Qk) ~ P (0) T a (q (0)) 



\l d (q,q k )-I(p(0),q(0))\ 



dL dL \ 

— {q,q) (!>! + — (q,q)<j) 1 dt + e q (q,q)j a (q k ) - p (0) T a (q (0)) 



h (dL 



o \dq 



(q,q) 



dL 



Q v (q (0) , q (0)) 0x (0) + e q (q, q) j a (q k ) - p (0) a (q (0)) 
Since q (t) solves the Euler-Lagrange equations, <fi\ (0) = 1 and <f»i (h) — 0, and q (0) = q k , 

T 



\l d (q,q k )-I(p(0),q(0)) \ = 



dL 

— (q (0) , q- (0))j a (q k ) + (e q (q, j)) 1 a (q k ) - p (0) 1 a (q k ) 
(p (0)) T a (q k ) + (e q (q, q)) T a(q k ) - (p (0)) T a (q k ) 

e q (q,q) a(qk) 



< \e q (q,q)\ \a(q k )\ 



<C a \e q {q,q)\ 

which yields the desired bound. 



□ 



Once again, if we assume that the quadrature error serves as a lower bound for the Sobolev error, combining 



the bounds from (3.3) and (3.4) yields 



\Id(q n ,qk) -I(p(t),q(t))\ < 2CC a (\\q n ~ q\\ w ^{[o,h]) + ll<?« ~ 9llz/»([o,ft]) + ||&> ~ «L«.([o,h])) ' 

This bound serves two purposes; the first is to establish a bound between the discrete Noether quantity and 
the Noether quantity computed on the local exact solution q. The second is to establish a bound between 
the discrete Noether quantity after one step and the Noether quantity computed on the initial data: 

\Id(q~n,qi) -i"O(0),g(0))| < 2CC a (\\q n - q\\ w i,i {[oM) + \\q„ - q~\\ L °°([o,h]) + M n ~ ^IL°°([o,ft])) ' 

since for (91,(72), q is the global exact flow of the Euler-Lagrange equations. 

The difference between these two bounds is subtle but important; by establishing a bound between the 
discrete Noether quantity and the Noether quantity associated with the initial conditions, on any step of the 
method we can bound the error between the discrete Noether quantity and the Noether quantity associated 
with the global exact flow. By establishing the bound between the discrete Noether quantity and the Noether 
quantity associated with q at any step, we can bound the error between the Noether quantity associated 
with the local exact flow q and the true Noether quantity conserved along the global exact flow: 

I J (p (t) , q (<)) - I (p (0) , q (0))| < I J (p (t) , q (t)) - I d (q n ,q k )\ + \I d (?„, q k ) - I (p (0) , q (0))| 

( 23 ) < ACC * (\\9n - q\\ w ^([o,h]) + Mn - g|lzao ([0jfc] ) + \\q n - ?|| £ oo ([0l h]) y 

for any to € [0, h] on any time step k. Because the local exact flow q is generated from boundary conditions 
(%, Qfe+i) which only approximate the boundary conditions of the true flow, there is no guarantee that the 
Noether quantity associated with q will be the same step to step, only that it will be conserved within 
each time step. However, because there is a bound between the Noether quantity associated with q and 
the discrete Noether quantity at every time step, the discrete Noether quantity and the Noether quantity 
associated with the exact flow, and because the Noether quantity is conserved point-wise along q on each 
time step, there exists a bound between the Noether quantity associated with each point of the local exact 
flow and the Noether quantity associated with the true solution. 

We finally arrive at the desired result, which is a theorem that bounds the error between the Noether 
quantity along the Galerkin curve and the true Noether quantity. It is significant because not only does it 
bound the error of the Noether quantity, but the bound is independent of the number of steps taken, and 
hence will not grow even for extremely long numerical integrations. 

Theorem 3.11. (Convergence of Conserved Noether Quantities) Define 

. dL , 

Pn = -qT {qn,q n ) ■ 



Under the assumptions of Lemmas (3.3 - 3.4), if the Noether map I (p,q) is Lipschitz continuous in both its 
arguments, then there exists a constant C v independent N , the number of method steps, such that: 

\I(p(0),q(0))- I(p n (t),q n (t))\ < C v (\\q n - q\\ w i, H [o, h]) + \\9n ~ q\\ L oo {[0th]) + \\q n - ?|| ioo([0i/l] )) ■ 
for any t G [0, Nh}. 

Proof. We begin by introducing the Noether quantity evaluated at t on the local exact flow, q: 

(24) I J (p (0) , q (0)) - I {pn (t) , q n (t))| < \I (pn (<) , <7« (*)) -Hp (*) , « (t))| 

+ |7(p(t),«(i))-/(p(0) > «(0))|. 
Considering the first term in p4| , let L4 be the Lipschitz constant for /(-,•). Then 
I J (p n (t) ,?„(*))-/ (p (t) ,q(t))\<U\ (p n (t) , q n (t)) - (p (t) ,q(t))\ 

<L 4 {\p n {t)~p{t)\ + \q n (t)-p(t)\) 



S(M*).9n(*))-|v (9 (*).?(*)) 



+ \q n (t)-q(t)\ 



< ^ (La |<7« (*) - g(t)| + (La + 1) l<7« (*) -?(*)!) 



(25) 



< L 4 (L 2 + 1) (||g„ - q\\ L ^ {[0Jl]) + \\q n 



iL°°([0,ti\) 



\q n 



lL°°([0,h]) 



+ 9n _ 9 



L°°([0 



,h])) 



The second term in (24) is exactly the bound given by (23) and thus combining (25) and (23) in (24) and 
defining C v = 4CC a + L (L 2 + 1), we have: 

\HPn (t),q n (t)) — I (p(t) ,q (t))\ < C v (\\q~n-q\\ w i.i([o,h]) 
which completes the result. □ 

The convergence and bounds of the Noether quantity evaluated on the Galerkin curve to that of the 
true solution is hampered by one issue. While Theorems (3.6) and (3.7) provide estimates for convergence 
in the Sobolev norm ||* II w rl « 1 ([o fc])' Theorem (3.11) requires estimates in the L°° norm. We can estab- 
lish a bound for \\q n (t) — q (i)llL°°([o /i])' but ^ ^ s muc h more difficult to establish a general estimate for 
Mn(t)-Ht)\\ L ~ {[0th] y 

Lemma 3.5. (Bound on L°° Norm from Sobolev Norm) For any t £ [0, h], the following bound holds: 

1 



\q (t)\ < max 



h' 



id thus 



W\\L<»([o,h]) 



< max 



1 



\1\\W 1 - 1 {10M]) 



Proof. This is a basic extension of the arguments from Lemma A.l. in Larsson and Thomee |11) . generalizing 
the lemma from the interval [0, 1] to an interval of arbitrary length, [0, h]. We note that for any t, s G [0, h], 
q(t) = q (s) + Jl q (u) dtt. Thus: 



\q(t)\<\q(s)\ 



\q(u)\du 



<!?(*)! + ll«llii([ ,fc]) 



Now, we integrate with respect to s: 

r-ll 



\q(t)\ds< / |g(«)|da+ / ||g|| L i ([ o,w) 
Jo JO 

h\q{t)\ < (||g|| £ i([ D| /,]) + '»ll9|li,i([o,fc])) 



, ds 



which yields the desired result. 



□ 



Under certain assumptions about the behavior of q n ~ q, it is possible to establish bounds on the point-wise 
error of q n from the Sobolev error \\q n — 9||tyi,i([o For example, if the length of time that the error is 
within a given fraction of the max error is proportional to the length of the interval [0,h], i.e. there exists 
C*i,C2 independent of h: i.e., 

m ({t||| (g„ (t) - Ut)) || > Ci || i (t)-f(«)|L})>«, 
where m is the Lebesque measure, then it can easily be seen that: 



\qn - ffllw^dOj/i]) 



> 



|9n (t) ~ q(t)\\ it > CiC 2 h \\q n - q\ 



L°°([0,Ji]) 



While we will not establish here that the q n converges in the L°° norm with the same rate that the Galerkin 
curve converges in the Sobolev norm, our numerical experiments will show that the Noether quantities tend 
to converge at the same rate as the Galerkin curve. 
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4. Numerical Experiments 



To support the results in this paper, as well as to investigate the efficiency and stability of spectral varia- 
tional integrators, several numerical experiments were conducted by applying spectral variational techniques 
to well-known variational problems. For each problem, the spectral variational integrator was constructed 
using a Lagrange interpolation polynomials at n Chebyshev points with a Gauss quadrature rule at 2n points. 
Convergence of both the one-step map and the Galerkin curves was measured using the £°° and L°° norms 
respectively, although we record them on the same axis using labeled L°° error in a slight abuse of notation. 
The experiments strongly support the results of this paper, and suggest topics for further investigation. 

4.1. Harmonic Oscillator. The first and simplest numerical experiment conducted was the harmonic os- 
cillator. Starting from the Lagrangian, 

L(q,q) = ^q 2 - \^ , 

where 9 £ R, the induced spectral variational discrete Euler-Lagrange equations are linear. Choosing the 
large time step h = 20 over 100 steps yields the expected geometric convergence, and attains very high 
accuracy, as can be seen in Figure [3j In addition, the max error of the energy converges geometrically, see 
Figure |4j and does not grow over the time of integration, see Figure [5] 

4.2. N-body Problems. We now turn our attention towards Kepler TV-body problems, which are both 
good benchmark problems and are interesting in their own right. The general form of the Lagrangian for 
these problems is 



N N i-1 



mirrij 



L(q,q) = lj2qfMq l + Gj2Y, 

i=l i=l 3=0 V" 

where qi € M. D is the center of mass for body i, G is a gravitational constant, and is a mass constant 
associated with the body described by q{. 

4.2.1. 2-Body Problem. The first experiment we will examine is the choice of parameters D — 2, mi = 
777*2 = L Centering the coordinate system at qi, we choose (72 (0) = (0.4,0), q^ (0) = (0,2), which has a 
known closed form solution which is a stable closed elliptical orbit with eccentricity 0.6. Knowing the closed 
form solution allows us to examine the rate of convergence to the true solution, and when solved with the 
large time step h = 2.0, over 100 steps, the error of the one-step map is O (0.56") with 77-refinement and 

with ^-refinement, as can be seen in Figure J^j and Figure j^J respectively. The numerical evidence 
suggests that our bound for the one-step map with /7-refinement is not sharp, as the convergence of the 
one-step map is always even. Interestingly, it is also possible to observe the different convergence rates of 
the one-step map and the Galerkin curves with n-refinement, as eventually the Galerkin curves have error 
approximately 0(0.74") while the one-step map has error approximately 0(0.56"), and V0.56 ~ 0.7483. 
However, it appears that the error from the one step map dominates until very high choices of 77, and thus 
it is difficult to observe the error of the Galerkin curves directly with /7-refinement, round off error becomes 
a problem before the error of the Galerkin curves does for smaller choices of n. 

The N-body Lagrangian is invariant under the action of SO(D), which yields the conserved Noether 
quantity of angular momentum. For the two body problem this is: 

L (q, q) = q x q v - q v q x 

where q = (q x , q y ). Numerical experiments show that the error of the angular momentum does not grow with 
the number of steps taken in the integration, Figure [13} but that the error is of the same order as the error of 
Galerkin curve with 77-refi nem ent in Figure [8j With /7-refinement, the angular momentum appears to have 



error O (h 2 +2) in Figure 11 This is interesting because the theoretical bound on the error of the Galerkin 
curves is O (h9), and the error of the Noether quantities is theoretically a factor C (h) times the error of the 
Galerkin curves, where C is the factor that arises in the proof of the convergence of the conserved Noether 
quantities. Numerical experiments suggest C is O (ft 2 ) for this problem, but that the Galerkin curves do 
converge at a rate of O (h?), which is consistent with of the Galerkin curve error estimate. While this 

24 



Convergence with N-Refinement 



10" 



10 



10 



t- 10 - 
o 



10"' 



10" - 



10 - 



9 


One Step Error 




Galerkin Curve Error 




e = 10 11 (0.21) N 



10" 



15 20 25 30 35 40 

Chebyshev Points Per Step (N) 



45 



Figure 3. Geometric convergence of the spectral variational integration of the harmonic 
oscillator problem, for 100 steps at step size h = 20.0. 

Energy Convergence with N-Refinement 
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Figure 4. Geometric convergence of the energy error of the spectral variational integration 
of the harmonic oscillator problem for 100 steps at step size h = 20.0. 
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Figure 5. Energy Stability of the spectral variational integration of the harmonics oscillator 
problem. This energy was computed for the integration using n = 14 for step size h = 20.0. 
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Figure 6. Geometric convergence of the Kepler 2-body problem with eccentricity 0.6 over 
100 steps of h = 2.0. Note that around n = 32, the error for the Galerkin curves becomes 
O (0.74™), while the error for the one-step map is always O (0.56™). 
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Energy Convergence with N-Refinement 
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Figure 7. Geometric convergence of the Energy Error of Kepler 2-body problem with 
eccentricity 0.6 over 100 steps of h = 2.0. Note that the error is O (0.74™), the same as it 
was for the Galerkin curves. 
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Figure 8. Geometric convergence of the angular momentum of the Kepler 2-body problem 
with eccentricity 0.6 over 100 steps of h = 2.0. Again, the error is of the same order as it 
was the Galerkin curves. 
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One Step Map Convergence with h-Refinement 
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Figure 9. Order Optimal convergence of the Kepler 2-body problem with eccentricity 0.6 
over 10 steps with h refinement. Note our bound is not sharp, as the error is 
where [•] is the ceiling function. 
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Figure 10. Convergence of the Kepler 2-body problem energy with eccentricity 0.6 over 
10 steps with h refinement. 
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Angular Momentum Convergence with h-Refinement 
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Figure 1 1 . Convergence of the angular momentum Kepler 2-body problem with eccentric- 
ity 0.6 over 100 steps of ft = 2.0. 
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Energy Behavior N=14 
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Figure 12. Stability of energy for Kepler 2-body problem. 
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Angular Momentum Behavior N=14 




FIGURE 13. Stability of angular momentum for Kepler 2-body problem. 



evidence is not conclusive, it is suggestive that the error analysis provides a plausible bound, 
analysis of the factor C would be an interesting direction for further investigation. 



A careful 



4.2.2. The Solar System. To illustrate the excellent stability proprieties of spectral variational integrators, 
we let D — 3, N — 10, and use the velocities, positions, and masses of the sun, 8 planet, and the dwarf 
planet Pluto on January 1, 2000 (as provided by the JPL Solar System ephemeris [23 ) as initial configuration 
parameters for the Kepler system. Taking 100 time steps of h = 100 days, the N = 25 spectral variational 
integrator produces a highly stable flow in Figure [14] It should be noted that orbits are closed, stable, and 
exhibit almost none of the "precession" effects that are characteristic of symplectic integrators, even though 
the time step is larger than the orbital period of Mercury. Additionally, considering just the outer solar 
system (Jupiter, Saturn, Uranus, Neptune, Pluto), and aggregating the inner solar system (Sun, Mercury, 
Venus, Earth, Mars) to a point mass, an N = 25 spectral variational integrator taking 100 time steps 
h = 1825 days (5 year steps) produces the orbital flow seen in Figure 15 Again, these are highly stable, 
precession free orbits. As can be clearly seen, the spectral variational integrator produces extremely stable 
flows, even for very large time steps. 

5. Conclusions and Future Work 

In this paper a new numerical method for variational problems was introduced, specifically a symplectic 
momentum-preserving integrator that exhibits geometric convergence to the true flow of a system under 
the appropriate conditions. These integrators were constructed under the general framework of Galerkin 
variational integrators, and made use of the global function paradigm common to many different spectral 
methods. 

Additionally, a general convergence theorem was established for Galerkin type variational integrators, 
establishing the important result that, under suitable hypotheses, Galerkin variational integrators will inherit 
the optimal order of convergence permitted by the underlying approximation space used in their construction. 
This result provides a powerful tool for both constructing and analyzing variational integrators, it provides 
a methodology for constructing methods of very high order of accuracy, and it also establishes order of 
convergence for methods that can be viewed as Galerkin variational integrators. It was shown that from 
the one step map, a continuous approximation to the solution of the Euler-Lagrange equations can be easily 
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Figure 14. Orbital diagram for the inner solar produced by the spectral variational inte- 
grator using all 8 planets, the sun, and Pluto with 100 time steps with h = 100 days. 




FIGURE 15. Orbital diagram for the inner solar produced by the spectral variational inte- 
grator using the 4 outer planets, Pluto, with the sun and 4 inner planets aggregated to a 
point with 100 time steps at h = 1825 days. 
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recovered over each time step. The error of these continuous approximations was shown to be related to 
the error of the one step map. Furthermore, the Noether quantities along this continuous approximation 
approximate the true Noether quantity up to a small error which does not grow with the number of steps 
taken. It was also shown that the error of the Noether quantities converges to zero with n or h refinement 
at a predictable rate. 

In addition to the convergence results, another interesting feature of spectral variational integrators is the 
construction of very high order methods that remain stable and accurate using time steps that are orders of 
magnitude larger than can be tolerated by traditional integrators. The trade off is that the computational 
effort required to compute each time step is also orders of magnitude larger than that of other methods, which 
are a major trade off in terms of the practicality of spectral variational integrators. However, a mitigating 
factor of this trade off is that the approach of solving a short sequence of large problems, as opposed to a 
large sequence of small problems, lends itself much better to parallelization and computational acceleration. 
The literature on methods for acceleration of the construction and solution of structured systems of linear 
and nonlinear problems for PDE problems is extensive, and it is likely that such methods could be applied 
to spectral variational integrators to greatly improve their computation cost. 

5.1. Future Work. Future directions for this work are numerous. Because of generality of the construction 
of Galerkin variational integrators, there exists many possible directions of further exploration. 

5.1.1. Lie Group Spectral Variational Integrators. Following the approach of Leok and Shingel [16] or Bou- 
Rabee and Marsden [1] , it is relatively straight forward to extend spectral variational integrators to Lie groups 
using natural charts. A systematic investigation of the resulting Lie group methods, including convergence 
and near conservation of Noether quantities, would be a natural extension of the work done here. 

5.1.2. Novel Variational Integrators. The power of the Galerkin variational framework is its high flexibility 
in the choice of approximation spaces and quadrature rules used to construct numerical methods. This 
flexibility allows for the construction of novel methods specifically tailored to certain applications. One 
immediate example is the use of periodic functions to construct methods for detecting choreographies in 
Kepler problems, which would be closely related to methods already used with great success to detect 
choreographies. Another interesting application would be the use of high order polynomials to develop 
integrators for high-order Lie group problems, such as the construction of Riemannian splines, which has 
a variety of applications in motion planning. Enriching traditional polynomial approximation spaces with 
highly oscillatory functions could be used to develop methods for problems with dynamics evolving on 
radically different time scales, which are also very challenging for traditional numerical methods. 

5.1.3. Multisymplectic Variational Integrators. Multisymplectic geometry (see Marsden et al. |20j ) has be- 
come an increasingly popular framework for extending much of the geometric theory from classical Lagrangian 
mechanics to Lagrangian PDEs. The foundations for a discrete theory have been laid, and there have been 
significant results achieved in geometric techniques for structured problems such as elasticity, fluid mechan- 
ics, non-linear wave equations, and computational electromagnetism. However, there is still significant work 
to be done in the areas of construction of numerical methods, analysis of discrete geometric structure, and 
especially error analysis. Galerkin type methods have become a standard method in classical numerical PDE 
methods, such as Finite-Element Methods, Spectral, and Pseudospectral methods. The variational Galerkin 
framework could provide a natural framework for extending these classical methods to structure preserving 
geometric methods for PDEs, and the analysis of such methods will rely on the notion of the boundary 
Lagrangian (see Vankerschaver et al. [55]), which is the PDE analogue of the exact discrete Lagrangian. 
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Appendix A. Proofs of Geometric Convergence of Spectral Variational Integrators 



As stated in {3.2 it can be shown that spectral variational integrators converge geometrically to the true 
flow associated with a Lagrangian under the appropriate conditions. The proof is similar to that of order 
optimality, and is offered below. 

However, before we offer a proof of the theorem, we must establish a result that extends Theorem |1.1| 
Specifically, we must show: 

Theorem A.l. (Extension of Theorem |1.1| to Geometric Convergence) Given a regular Lagrangian L and 
corresponding Hamiltonian H, the following are equivalent for a discrete Lagrangian Ld (qo,qi,n): 

(1) there exist a positive constant K , where K < 1, such that the discrete Hamiltonian map for Ld (qo, qh, n) 
has error O (K n ), 

(2) there exists a positive constant K, where K < 1, such that the discrete Legendre transforms of 
L d (qo, qh, n) have error O (K n ), 

(3) there exists a positive constant K, where K < 1, such that Ld (qo, Qhi n ) approximates the exact 
discrete Lagrangian L^ (qo, qh, h) with error O (K n ). 

The proof of this theorem is a simple modification of the proof of Theorem and is included here for 
completeness. For details, the interested reader is referred to |19j . 

Proof. Since we are assuming that the time step h is being held constant, we will suppress it as an argument 
to the exact discrete Lagrangian, writing (qo, qh) for L^ (qo, qh,h). First, we will assume that Ld (qo, qh, n) 
approximates Ld (qo, qh) with error O (K n ) and show this implies the discrete Legendre transforms have error 
O (K n ). By assumption, if Ld (qo,qh,n) has error O (K n ), there exists a function which is smooth in its first 
two arguments e v : Q x Q x N — >• R such that: 

Ld (qo,qh,n) = (qo,Qh) + K n e v (q ,qh,n) , 

with \e v (qo,qh,n)\ < C v on U v . Taking derivatives with respect to the first argument yields: 

F~£3 (g , qh) = ^~L% (qo, q h ) + K n D ie . v (q , q h , n) , 
and with respect to the second yields: 

¥+L n d (qo, q h ) = F+if (go, qh) + K' l D 2 e v (q , q h ,n). 

Since e v is smooth and bounded over the closed set U, so are Die v and D 2 e v , yielding that the discrete 
Legendre transforms have error O (K n ). Now, to show that if the discrete Legendre transforms have error 
O (K n ), the discrete Lagrangian has error O (K n ), we write: 

e-v (qo,Qh,n) = — [L d (q a ,qh,n) - if (q ,qh)] , 

L>ie v (q , q h , n) = ^ [F~ L d (q , q h , n) - F^Lf (qo,q h )] , 

L> 2 e v (q ,qh,n) = — - [¥ + L d (qo,qh,n) - F+if (qo,qh)] ■ 

Since D\e v and D 2 e v are smooth and bounded on a bounded set, this implies there exists a function d(n) 
such that 

\\e v (q(0),q(h),n)-d(n)\\ < C v , 

for some constant C v . This shows that the discrete Lagrangian is equivalent to a discrete Lagrangian with 
error O (K n ). We note that the equivalence is a consequence of the fact that one can add a function of h 
or n to any discrete Lagrangian and the resulting discrete Euler-Lagrange equations and discrete Legendre 
Transforms are unchanged, hence the function d(n). 

To show the equivalence of the discrete Hamiltonian map having error O (K n ) and the discrete Legendre 
transforms having error O (K n ), we recall expressions for the discrete Hamiltonian map for the discrete 
Lagrangian Ld and exact discrete Lagrangian L?' 

F Ld =¥+L d o (¥-L d y\ 
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F l k = F+Lf o (F-Lf ) \ 

Now, we make use of the following consequence of the implicit function theorem: If we have smooth functions 
gi,g 2 and the sequences of functions {/i„}^Li, {f2 n }n=v { e i«}^=i and { e 2„}^Li sucn that 

/i„ (z) = ff i (as) + X"e u (a) , 

/ 2 „ (z) = 92 (x) + K n e 2n (x) , 
where sup{||ei n ||}^_ 1 < C\ and sup{||e2„||}^ 1 < C 2 on compact sets, then 

(26) f 2n (fi n (x)) = .92 (3i (x)) + K n e 12n (x) 

(27) fi n 1 (v)=9i 1 (y)+K n e ln (y) 

for some sequences of functions {ei2„}^L 1 , {ei„}^Li where sup{||ei2„||}^L 1 < C\ and sup {||ei n \\}™ = i < C 2 
on compact sets. 

It follows from (26) and (27) that if the discrete Legendre transforms have error 0(K n ), the discrete 
Hamiltonian map does as well. Finally, if we have a discrete Hamiltonian map has error O (K n ), we use the 
identity 

(¥~L d ) 1 (q ,po) = {qo,KQ o F Ld {q ,p )) 

where ttq is the projection map ttq : (q,p) —> q and (27) to see that ¥~L d is has error O (K n ), and the 
identity: 



¥+L d = F Ld o¥-L d , 



along with (26) to establish that F + L d also has error O (K n ), which completes the proof. 



□ 



This simple extension is a critical tool for establishing the geometric convergence of spectral variational 
integrators, and leads to the following theorem concerning the accuracy of spectral variational integrators. 

Theorem A. 2. (Geometric Convergence of Spectral Variational Integrators) Given an interval [0, h] and a 
Lagrangian L : TQ — > R, let q be the exact solution to the Euler- Lagrange equations, and q n be the stationary 
point of the spectral variational discrete action 



Ll(q ,q h ,n)= ^ext^ S d ({«}^) 

<ln(0)=qo,q n (h)=q h 



ext h V* b n , L (q n (c„ h) , q n (c n h) ) 

qn(0)=q ,q n (h)=q h 



If: 



(1) there exists constants Ca,Ka, Ka < 1, independent ofn, such that, for each n, there exists a curve 
q n E M™ ([0,h],Q) such that, 



(9,0) - Un,Qn) 



<C A K n A 



A- 



(2) there exists a closed and bounded neighborhood U C TQ such that (q (t) ,q(t)) € U and (^q n (f) , q n (t) 
U for all t and n, and all partial derivatives of L are continuous on U , 

(3) for the sequence of quadrature rules Q n (/) = /i^"^ b nj f (c nj /i) ~ J Q f (t) dt there exists constants 
C g , Kg, Kg < 1, independent ofn such that: 



L (q n (t) ,q n (t)) dt - h^b nj L (q n (c nj h) ,q n (c nj h)) 

3 = 1 

for anyq n eM n {[0,h},Q), 
(4) and the stationary points q, q n minimize their respective actions, 



< CgK-. 



then 
(28) 



\Lf (q ,qi)-L s d (q ,q u n)\<C s K2 

for some constants C S1 K S , K s < 1, independent of n, and hence the discrete Hamiltonian flow map has 
error O (K™). 
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Proof. As before, we rewrite both the exact discrete Lagrangian and the spectral discrete Lagrangian: 

rh 

L (q (t) , $ (*)) dt - G n (L (q n (t) , (t))) 



l-kf (<?o,9i) - (go, Si, «) I 



h) An (Cnjh)) 



j = l 

L (q, q)dt-h^2 

Vrij ^ [Qn, Qn 

3 = 1 

with suppression of the t argument. We introduce the action evaluated on the curve q n : 
(29) 



L(q,q) dt-h^2 

&nj ^ [Qn, Qn 

J'=l 



L(q,q)dt- / L 



dt 



(30) 



(31) 



< 



L [qn , gn j dt - ft ^ 6 nj L (g n , <?„) 
i=i 

I (r/. </)<!/ - / L dt 

/ L ( q n , g n J dt - ft ^ 6„ 3 L (q n , g„ 
■/o 7=1 



Considering the first term in (31): 



J L(q,q)dt- J L(q n ,qn^dt 



< 



(Q,Q) - L (q n ,q"n) 
(Q,Q) - L (q n An) 



dt 



dt. 



By assumption, all partials of L are continuous on U, and since U is closed and bounded, this implies L 
is Lipschitz on U , so let L a denote the Lipschitz constant. Since, again by assumption, (q, q) € U and 

Qn, Qn € C/, we can obtain: 



L 



(?) Q)~ L (q n ,qn^ 



dt < / L Q 



(«,?) - (<Zn,<7n) 



dt 



< / L a C A K\dt 
Jo 

=hL a CAK^. 



Hence, 
(32) 







L(q,q)dt- L(q n ,qn)dt 



< hL a C A K r , 



Next, considering the second term in (31 ) 

i-h 



/ L (q„An) dt - S^hb^L (q„An) 

J ° V 7 , = 1 



36 



since q n minimizes its action 
(33) 



b n?i L (grugn) < b n .L ( gn,gnj < / L ( q n , q„) dt + CgKg 

.7 = 1 3=1 J ° 



where the inequalities follow from the assumptions on the order of the quadrature rule and ( 32 ). Furthermore, 

hJ2b nj L(q n ,t n ) > L (q n , q n ) dt - C g K g > / L (q,q) dt - C g K n g 
j=1 Jo Jo 



(34) 



' J L(q n , kn) dt - hL a C A K\ - C g K n g , 



where the inequalities follow from (32), the order of the sequence of quadrature rules, and the assumption 
that q minimizes its action. Putting (33) and (34) together, we can conclude: 



(35) 



/ L [q n ,g„j dt-hy] b n . L (q n , q„ 

^° 7 = 1 



< (hL a C A + C g ) K, 



— n 

9 1 xl s ■ 



where K s — max {K A , K g ). Now, combining the bounds (32) and (35) in (31), we can conclude 

\L E d (q Q , qi )-L s d (q ,qi,n)\ < (2hL a C A + C g ) K~ n 
which, combined with Theorem |A.1| establishes the rate of convergence. 



□ 
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